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Abstract 

In this paper, a numerical method for the solution of a strongly coupled reaction-diffusion system, with suitable initial and 
Neumann boundary conditions, by using cubic B-spline collocation scheme on a uniform grid is presented. The scheme is 
based on the usual finite difference scheme to discretize the time derivative while cubic B-spline is used as an interpolation 
function in the space dimension. The scheme is shown to be unconditionally stable using the von Neumann method. The 
accuracy of the proposed scheme is demonstrated by applying it on a test problem. The performance of this scheme is 
shown by computing L x and Li error norms for different time levels. The numerical results are found to be in good 
agreement with known exact solutions. 
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Introduction 

This study is concerned with the numerical solution of strongly 
coupled reaction-diffusion system using cubic B-spline. Reaction- 
diffusion system arises in the study of biology, chemistry and 
population dynamics. It can be used to describe a mathematical 
model of a class of chemical exchange reaction that arises in the 
transport of ground water in an aquifer [1]. The mathematical 
formulation of the problem is: 



8 /u 

Ttlv 



a 0 

P y 



+ 



f(u,v) 



xe[a,b], te[0,T] 



with the initial conditions 



u(x,Q) = uo(x), v(x,0) = vo(x), 
and the boundary conditions 

u x (0,t) = ui(i), u x {\,t) = p l {t), 

v. Y (0,r) = s( 2 (0, v. v (U) = &(0, 



(1) 



(2) 



(3) 



where u = u(x,t) and v = v(x,f) are the concentrations of the two 
substances in the interaction and the constants a,/8,y are such that 
a > 0,P > 0, y 0. The following consistency conditions hold 



oc 1 (0) = M' 0 (0), p l (0) = u' 0 (\), 
oc 2 (0) = v'o(0), fe(0) = v' 0 (l). 



(4) 



The global solutions of such a system of equations have attracted 
the attention of several researchers [2-6]. Researchers have also 
investigated the existence, uniqueness and boundedness of the 
global solution in bounded and unbounded region [3,6]. Cao and 
Sun [1] derived a finite difference scheme by the method of 
reduction of order for the numerical solution of strongly coupled 
reaction-diffusion system with Neumann boundary values condi- 
tions. They proved the solvability and convergence by using the 
energy method. Several researchers focused on analytical solutions 
of nonlinear equations by using approximate analytical methods. 
Examples include Ghoreishi et al [7] who obtained the analytical 
solution for a strongly coupled reaction-diffusion system by using 
the Homotopy Analysis Method. The solution of the system was 
calculated in the form of an infinite series with easily computed 
components. This method cannot always guarantee the conver- 
gence of approximate series [8,9] and the method depends on 
choosing the proper linear operator and initial guesses. 

The study of spline functions is a key element in computer aided 
geometric design [10] and also several other applications. It has 
also attracted attention in the literature for the numerical solution 
of linear and non-linear system of second-order boundary value 
problems [11,12] that arise in science and engineering. Some 
researchers have considered spline collocation method for 
diffusion problems [15-19]. Advection-diffusion equation arises 
frequently in the study of mass, heat, energy and vorticity transfer 
in engineering. Bickley [13] introduced the idea of using a chain of 
low-order approximation (cubic splines) rather than a global high- 
order approximation to obtain better accuracy for a linear 
ordinary differential equation. Fyfe [14] used the method 
proposed by Bickley [13] and conducted an error analysis. It 
was concluded that the spline method is better than the usual finite 
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difference scheme because it has the flexibility to obtain the 
solution at any point in the domain with greater accuracy. The 
numerical solution of some partial differential equations can be 
obtained using B-spline functions of various degrees which can 
provide simple algorithms. As an example, a combination of finite 
difference approach and cubic B-spline method was used to solve 
the Burgers' equation [15], heat and wave equation [16,17], 
advection-diffusion equation [18] and coupled viscous Burgers' 
equation [19]. 

Sahin [24] presented the B-spline methods in several degrees for 
the solution of following non linear reaction-diffusion system 



8U_ 
8t 



dV 
lit' 



d 2 U 
Jx 2 



+ b ] U + c ] V + d l U 2 V + e l UV + m 1 UV 2 + n u 



c- 1 V + d 1 U 2 V + e 1 UV + m 2 UV 2 +n 1 . 



The finite element method was employed for the solution of 
reaction-diffusion systems and the time discretization of the system 
was achieved by the using Crank-Nicolson formulae. The system 
was considered only reaction-diffusion but the problem we 
propose to investigate is a strongly reaction-diffusion system with 
an extra term U xx in the second equation of system (1). 

In this paper we aim to apply the combination of finite 
difference approach and cubic B-spline method to solve the system 
(1). Some researchers have utilized the B-spline collocation 
methods to solve systems of differential equations but so far as 
we are aware not the system (1). A usual finite difference scheme is 
used to discretize the time derivative. Cubic B-spline is applied as 
an interpolation function in the space dimension. The uncondi- 
tional stability property of the method is proved by von Neumann 
method. The feasibility of the method is shown by a test problem 
and the approximated solutions are found to be in good agreement 
with the exact solution. 

This paper is structured as follows: Firstly, we discuss a 
numerical method incorporating a finite difference approach with 
cubic B-spline. Secondly, the von Neumann approach is used to 
prove the stability of method and compare the numerical solution 
with exact solution of system (1). Thirdly, a test problem is 
considered to show the feasibility of the proposed method. Finally, 
the conclusion of this study is given. 

Description of Cubic B-Spline Collocation Method 

In this section, we discuss the cubic B-spline collocation method 
for solving numerically the strongly coupled reaction-diffusion 
system (1). The solution domain a<x<b is equally divided by 
knots Xi into n subintervals £=0,1,2,. ..,« — 1 where 

a = Xq <X\ < ... <x„ = b. Our approach for strongly coupled 
reaction-diffusion system using collocation method with cubic B- 
spline is to seek an approximate solution as [20] 



u i (x,t)=j2 , ;: l _ 3 c,{t)B 3 , i ( x ), 



(5) 



where C,(/) and D t (t) are (time dependent) quantities which are to 
be determined for the approximated solutions Ui(x,t) and Vi(xj) 
to the exact solutions u(x,t) and v(x,t) respectively, at the point 
(xi,tj) and B 3i {x) are cubic B-spline basis functions which are 
defined by [21] 



1 



P + V?(x -X,+ i) + 3h(x-x i+l f - 3(x-X i+ i ) 3 
A 3 + 3A 2 (x i+ 3 - x) + 3h(x,+ 3 - xf - 3(x i+3 - xf 

(x i+ 4-xf 

0 



XE[Xi,X i+i \ 

xe[x i+ uXi+z] 
xe[xt+z,x i+ 3,] 
xe[Xi+3,Xt+t] 
otherwise 



(6) 



where h = (b — a)/n. The approximations U\ and V\ at the point 
(Xi,tj) over subinterval + can be defined as 



vi- 



(7) 



where i = 0,1,2,. ..,». So as to obtain the approximations to the 
solutions, the values of 2?3 ; ;(x) and its derivatives at nodal points 
are required and these derivatives are tabulated in Table 1. 

Using approximate functions (6) and (7), the values at the knots 
of Uj and Vj and their derivatives up to second order are 
determined in the terms of time parameters Ci and Z)^ as 



1 2 1 

tjj = _c' +-C J +~C J 

1 



i 



1 2 i 

1 2 1 

<T Y = — D 1 --D 1 +—D' 



(8) 



h 2 



The approximations of the solutions of the system (1) at tj' +l time 
level can be given by as [22] 



(E/,X: = 6$ +1 +(l-0)/i 



(9) 



where />', 7(<\ x ):-//. = t/^+KM+Si and 
fj = a 3 C// + tt. 4 V J i , g' i = oi 5 U' j +a b Vj, where a,-,i'= 1,2,3,4 are con- 
stants and the subscripts j and j + 1 are successive time levels, 

7 = 0,1,2, Discretizing the time derivatives in the usual finite 

difference way and rearranging the equations, we obtain 



7+1 
i 



■ Ui + (l-9)At I / i , 
:Vl + (l-0)At<ji, 



(10) 



Table 1. Values of Bxi(x) and its derivatives. 









Xm 










0 


1/6 


2/3 


1/6 


0 


B'i 


0 


l/2/l 


0 


-1/2h 


0 


B"i 


0 


1/h 2 


-21 h 2 


1/h 2 


0 
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Table 2. Some numerical solution of U(x,f) at t=\.0. 





nix 


0.125 


0.25 


0.375 


0.625 


0.75 


0.875 


16 


0.09166515 


0.18393918 


0.27621321 


0.27621321 


0.18393918 


0.09166515 


32 


0.06452148 


0.18393899 


0.30335651 


0.30335651 


0.18393899 


0.06452148 


64 


0.05661 594 


0.18393653 


0.31125712 


0.31125715 


0.18393658 


0.05661601 


128 


0.05457118 


0.18394479 


0.31331839 


0.31331836 


0.18394472 


0.05457110 


256 


0.05407447 


0.18393318 


0.31385834 


0.31385863 


0.18396692 


0.05407518 


512 


0.05392663 


0.18393549 


0.31405060 


0.31405084 


0.18394405 


0.05392761 


1024 


0.05386625 


0.18393840 


0.31400978 


0.31400981 


0.18394110 


0.05388334 


"exactly 1) 


0.05387469 


0.18393972 


0.31400474 


0.31400474 


0.18393972 


0.05387469 
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where At is the time step. It is noted that the system becomes an 
explicit scheme when 8 = 0, a fully implicit scheme when 0 = 1, 
and a Crank-Nicolson scheme when 8 = 0.5 [15] with time 
stepping process being half explicit and half implicit. In this paper, 
we use the Crank-Nicolson approach. Hence, (10) becomes 



t/ +1 -0.5A/^ 



-0.5Af^ 



■ U{- 
V{- 



•0.5A^, 
0.5A*^. 



(11) 



From equations (11) and (12), the system becomes a matrix system 
of dimension 2(w + 3) x 2(« + 3) which is a bi-tridiagonal system 
that can be solved by the Thomas Algorithm [23]. From equation 
(10) and (12), the system can be written in the matrix vector form 
as: 



ME' + I =NE' + b, 



(13) 



The system thus obtained on simplifying (11) after using (8) consists 
of 2(«+l) linear equations in 2(m + 3) unknowns 

o '=(d+ 3 l ,d+ 2 l .C1+ 1 ,...,d„ + _ i) , \y=(p>^ ,d<'+ 2 1 ..,D> n +ty 

at the time level t = tj+\. So as to obtain a unique solution to the 
resulting system, four additional linear equations are required. 
Thus, equation (7) is applied to the boundary conditions (3) to 
obtain 



where 



£'=[c': 3 ,c: 2 ,c^ 1 ,...,c^ 1 ,iy:3,iy: 2 ,/y: 1 ,...,iy„_ 1 ] r , 



(«74 +1 =ai(f,- + i), 

M +1 =ft(*,+i), 

(Vj 0 +i =« 2 (tj +l ), 

(v x y„ +1 =fs 2 (t J+1 ), 



(12) 



b = {oc l (t j+l )&o,...Mtj+i)Mt J+ i)m---Mtj+i)} J=o,h2, 

and M is an 2(« + 3) x 2(« + 3) dimensional matrix given by 



Table 3. Some numerical solution of V{x,t) at t= 1.0. 



nix 


0.125 


0.25 


0.375 


0.625 


0.75 


0.875 


16 


0.27700610 


0.18393970 


0.09087339 


0.09087339 


0.18393970 


0.27700601 


32 


0.30361313 


0.18393970 


0.06426627 


0.06426627 


0.18393970 


0.30361313 


64 


0.31132848 


0.18393971 


0.05665094 


0.05665093 


0.18393969 


0.31132845 


128 


0.31333062 


0.18393969 


0.05454876 


0.05454877 


0.18393971 


0.31333065 


256 


0.31383601 


0.18393980 


0.05404358 


0.05404347 


0.18393960 


0.31383575 


512 


0.31396262 


0.18393980 


0.05391696 


0.05391688 


0.18393963 


0.31396241 


1024 


0.31399647 


0.18394147 


0.05388465 


0.05388429 


0.18393793 


0.31399625 




0.31400474 


0.18393972 


0.05387469 


0.05387469 


0.18393972 


0.31400474 
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M-- 



A 2 
A, 



A 5 



N-- 



A 6 



po 0 -p 0 0 ... 0 0 
pi q x pi 0 ... 0 0 
0 p\ q x p\ 0 ... 0 



A 2 = 



0 0 

0 0 

0 0 0 0 

pi qi pi 0 

0 p 2 q 2 Pi 



0 



pi q\ pi 
Po 0 -p 0 

0 0 " 

0 0 

... 0 



0 0 
0 0 



pi qi pi 

0 0 0 



Ay 



0 0 0 0 ... 0 0 
Pi q 5 p 5 0 ... 0 0 
0 p 5 q 5 p 5 0 ... 0 



A b = 



0 0.. 

0 0.. 

0 0 0 0 

Pb qb Pb 0 

o Pb qb Pb 



ps qs ps 



0 



0 0 
0 0 



Pb qb Pb 
ooo 



A 3 -- 



Po 


0 


~Po 


0 


Pi 


93 


Pt, 


0 


0 


Pi 


q-i 


P3 



A 4 = 



0 0 

0 0 

0 0 0 0 

Pa q<i P4 0 

0 p 4 q 4 p 4 



Pi qi Pi 
Po 0 -p 0 

... 0 0 

... 0 0 

0 ... 0 



0 0 
0 0 



Pa 94 P4 
0 0 0 



Ar- 



0 0 0 0 ... 0 0 

pi q-] pi 0 ... 0 0 
0 pi qi pi 0 ... 0 



A s = 



0 0.. 

0 0.. 

0 0 0 0 

Ps q% Pt, 0 

0 Ps qi Pi 



pi qi pi 

0 0 0 _ 

... 0 

... 0 



0 



0 0 
0 0 



Pi qs Pi 

0 0 0 
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Table 4. The absolute errors at different grid points of 
numerical solution of U(x,t) at t=l.O. 


Table 5. The absolute errors at different grid points of 
numerical solution of V(x,t) at f = 1.0. 






n/x 


0.125 0.25 0.375 0.625 0.75 0.875 


n/x 


0.125 0.25 0.375 0.625 0.75 0.875 


16 


3.7790E-02 5.3228E-07 3.7791 E-02 3.7791 E-02 5.3354E-07 3.7790E-02 


16 


3.6998E-02 1.5559E-08 3.6998E-02 3.6998E-02 1.5097E-08 3.6998E-02 


32 


1.0646E-02 7.2246E-07 1.0648E-02 1.0648E-02 7.2208E-07 1.0646E-02 


32 


1.0391 E-02 1.5257E-08 1.0391 E-02 1.0391 E-02 1.5398E-08 1.0391 E-02 


64 


2.7412E-03 3.1889E-06 2.7476E-03 2.7476E-03 3.1346E-06 2.7413E-03 


64 


2.6762E-03 5.3785E-09 2.6762E-03 2.6762E-03 2.5277E-08 2.6762E-03 


128 


6.9649E-04 5.0748E-06 6.8634E-04 6.8638E-04 5.0093E-06 6.9640E-04 


128 


6.7411E-04 2.7316E-08 6.7406E-04 6.7407E-04 3.3408E-09 6.7408E-04 


256 


1.9977E-04 6.5450E-06 1.4640E-04 1.4610E-04 2.7209E-05 2.0048E-04 


256 


1.6873E-04 8.3759E-08 1.6888E-04 1.6877E-04 1.1441 E-07 1.6898E-04 


512 


5.1944E-05 4.2320E-06 4.5862E-05 4.6101E-05 4.3318E-06 5.2920E-05 


512 


4.2123E-05 8.0580E-08 4.2271 E-05 4.2184E-05 8.0882E-08 4.2334E-05 


1024 


8.4430E-06 1.3219E-06 5.0421E-06 5.0741E-06 1.381 1E-06 8.6510E-06 


1024 


8.2681 E-06 1.7575E-06 9.9601 E-06 9.5943E-06 1.7881 E-06 8.4926E-06 


doi:1 0.1 371 /journal.pone.0083265.t004 
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2yAt6 t 2(l-<x 6 At0) 

Also the entries in sub-matrices Aj,i= 1,2,3,. ..,8 have the ?3 = 

following form 



h 2 



a.At6 (l-a 3 A*0) (BAt6 ct 5 At6\ 



laiAtO 2(l-a 3 A*0) 
~1P~ + 3 ' 



P2- 



?2 = 



P3- 



yAt6 (l-a 6 A«0) 

~h r+ 6 ' 



?4 = 



/2/JAr0 2oc s At6\ 



\ h 2 



I' 



P5- 



(l + a 3 A«(l-0)) aAf(l- 



h 2 



?5 = 



2(l + a 3 Af(l-0)) 2aA/(l-0) 
3 h 2 ' 



P6- 



a 4 Af(l-0) 
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Figure 2. The absolute error curves of numerical solution V(x,t) at /-l.O. 
doi:1 0.1 371 /journal.pone.0083265.g002 



96 = 



2a 4 A<l-0) 



(l + a 6 A?(l-0)) yAf(l- 

P7 = 6 + 



97 = 



2(l + a 6 A/(l-0)) 2yAi(l-0) 
3 F ' 



a particular time level can be calculated repeatedly by solving the 
recurrence relation [19]. 

C° and D° can be obtained from the initial condition and 
boundary values of the derivatives of the initial condition as follows 
[15]: 



(U?) x = u' 0 ( Xl ),i = 0, 
Uf = u 0 (xi),i=0,l,2,...,n, 
( u ?) x = u 'o(x i ),i = n, 



(14) 



ftAt(l-9) Ks Af(l- 



h 2 



and similarly for approximate solution V\ 



+ 1 



?8 = 



2jSA*(l-0) 2a 5 A/(l- 



/? 2 



Initial state 

After the initial vectors C° and D° have been computed from 
the initial conditions, the approximate solutions U{ + and Vj +1 at 



(Vf) x = v' 0 ( Xi ),i=0, 

V? = vo(xi),i=a,l,2,...,n, 

(Vf) x = v'o(xd,i=n. 



(15) 



Thus the equations (14) and (15) yield a 2(w + 3) x 2(« + 3) matrix 
system, of the form 

AEi +1 =d 



Table 6. The absolute errors of numerical solution [1] of 
u(x,t) at t= 1.0. 



nix 


0.125 


0.375 


0.625 


0.875 


16 


1.3906E-03 


9.7932E-04 


9.7932E-04 


1 .3906E-03 


32 


3.471 9E-04 


2.4432E-04 


2.4432E-04 


3.471 9E-04 


64 


8.6706E-05 


6.0986E-05 


6.0986E-05 


8.6706E-05 


128 


2.1670E-05 


1 .5240E-05 


1 .5240E-05 


2.1670E-05 


256 


5.41 73E-06 


3.8098E-06 


3.8098E-06 


5.41 73E-06 


512 


1.3543E-06 


9.5242E-07 


9.5242E-07 


1.3543E-06 
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Table 7. The absolute errors of numerical solution [1] of 



v(x 


,t) at /= 1.0. 










nix 


0.125 


0.375 


0.625 


0.875 


16 


1.1016E-03 


1.5129E-03 


1.5129E-03 


1.1016E-03 


32 


2.6765E-04 


3.7052E-04 


3.7052E-04 


2.6765E-04 


64 


6.6834E-05 


9.2554E-05 


9.2554E-05 


6.6834E-05 


128 


1.6703E-05 


2.3133E-05 


2.3133E-05 


1 .6703E-05 


256 


4.1756E-06 


5.7831 E-06 


5.7831 E-06 


4.1756E-06 


512 


1.0438E-06 


1.4457E-06 


1.4457E-06 


1 .0438E-06 
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Table 8. Errors at different time-levels for U(x,t) with 
A/ = 0.001. 





n=100 




n=300 




/7=500 




f 


L2 


U 




L_ 




U 


0.01 


5.90E-05 


6.95E-05 


6.56E-06 


7.67E-06 


2.36E-06 


2.76E-06 


0.1 


4.24E-04 


5.47E-04 


4.71 E-05 


6.04E-05 


1 .69E-05 


2.17E-05 


0.5 


1.33E-03 


2.56E-03 


1 .48E-04 


2.83E-04 


5.34E-05 


1.01 E-05 


1.0 


1 .60E-03 


5.06E-03 


1.81E-03 


5.61 E-03 


8.10E-05 


2.15E-05 



doi:1 0.1 371 /journal.pone.0083265.t008 

where 



Pa 


0 


~Pa 


0 






0 


0 


0 


0 


0 






0 ' 


P 


? 


P 


0 






0 


0 


0 


(1 


0 






() 


0 


P 


'/ 


p 






0 


0 


0 


(1 


0 






() 


0 


0 






p 


« 


p 


0 


0 






0 


0 


0 


0 


0 






Pa 


0 


-pa 


0 


0 






() 


() 


0 


0 


() 


0 


0 






0 


Pa 


0 


-Pa 


0 






() 


0 


0 


() 


0 






0 


P 


9 


!' 


0 






0 


0 


0 


0 


0 






0 


0 


/' 


9 


p 






0 


0 


0 






0 


0 


0 


0 


0 






/» 


'/ 


p 


0 


0 






0 


0 


0 


0 


0 






/'II 


0 


-Pa. 



The solution of above system can be found by Thomas 
Algorithm. 

Stability of the Scheme 

In this section, von Neumann stability method is applied for 
investigating the stability of the proposed scheme. This approach 
has been used by many researchers [15-19]. Substituting the 
approximate solution U and V, their derivatives at the knots and 
f(u,v)=(fif 2 ), g(u,v)=(gugz) into equation (10) yields a 
difference equation with variables C m and D,„ given by: 



«5(Ci ! _ 3 + C4_ 1 )+« 6 C{,_ 2 +a 7 (iy m _3 
"»(<£-3 + d+\ ) + «io C<:\ + an 



-r/'+l 



V+ 1 



where 



-^_i)+^ m _ 2 , 

(16) 



"2 = 



(1- Are/1 ) Arfe 
6 A 2- ' 



2(1- A/0/! ) 2Ar0a 

3 



Cl3 z 



At8f 2 



E j+\, 
7 = 0, 



=[c /+ 3 i ,c / : + 2 i ,c: + I 1 , 



W+l 7j/+l ,y+l jy + l ni+ l ] T 



= [Ko(x 0 ),i/o(^o),«o(-«l ),..., 

Uo(x„),u' 0 (x n ),V 0 (x 0 ),v 0 (x 0 ),v 0 (x l ),..., 
Vo(x„),v' 0 (x„)] T , 



A 1 1 2 

and„ 0 =^=^=3. 



«4 = 



2At6f 2 
~3~' 



a$ = 



(l + A?(l-0)/i) A?(l-0)a 
_ 6 + ■ 



«6 = 



2(l + A/(l-0)/i) 2A*(l-0)a 

^3 # ' 



a 7 = 



Table 9. Errors at different time-levels for V(x,t) with 
A/ = 0.001. 





n=100 




n = 300 




n = 500 




t 






L 2 


L_ 






0.01 


1 .05E-05 


1 .22E-05 


1.16E-06 


1.36E-06 


4.21 E-07 


4.90E-07 


0.1 


3.25E-04 


4.14E-04 


3.62E-05 


4.61 E-05 


1.30E-05 


1 .66E-05 


0.5 


1 .26E-03 


2.40E-03 


1 .40E-04 


2.68E-04 


5.07E-05 


9.65 E-05 


1.0 


1.55E-03 


4.87E-03 


1 .74E-04 


5.45E-04 


6.39E-05 


1 .96E-05 



doi:1 0.1 371 /journal.pone.0083265.t009 



"8 = 



2A/(l-0)/ 2 



/A?6gi Atda 



2A/0a 2ArO gi 
~h 2 3 
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V.(>Lt) 




Figure 3. A comparison between numerical U(x,f) and exact solutions at different time levels. 

doi:1 0.1 371 /journal.pone.0083265.g003 



(1-At6g 2 ) At8p (l-0)A/gi (l-0)Afa 
= 6 W au= 6 + W ' 




0.: 0.4 0.6 0.8 1.0 

Figure 4. A comparison between numerical V(x,t) and exact solutions at different time levels. 

doi:1 0.1 371 /journal.pone.0083265.g004 
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Table 10. Maximum error, ratio and order of convergence of 
the proposed scheme for U(x,t). 



doi:1 0.1 371 /journal.pone.0083265.t01 0 



(l + (l-fl)A/g 2 ) (\-8)Atp 
6 h 2 



f=0.1 






Order of 


f=0.5 




Order of 


n 


U 


Ratio 


Conver. 


to 


Ratio 


Conver. 


16 


1 .64E-02 






4.8E-02 






32 


4.13E-03 


3.9623 


1.9863 


1 .27E-02 


3.7704 


1.9147 


64 


1 .04E-03 


3.9906 


1.9966 


3.24E-03 


3.9408 


1 .9785 


128 


2.59E-04 


3.9975 


1.9991 


8.13E-04 


3.9850 


1 .9945 


256 


6.48E-05 


3.9990 


1.9996 


2.03E-04 


3.9959 


1 .9985 


512 


1.62E-05 


3.9984 


1 .9994 


5.09E-05 


3.9983 


1 .9994 


1024 


4.06 E-06 


3.9973 


1 .9979 


1 .27E-05 


3.9868 


1 .9952 



x 2 



A-m.Afi-m 



(cosi/f/i+2) + 



IttfBAd 
h 2 



' 1 — cos \j/h) 



X, 



B+At(\-8)(A gl +B g2 j 



(cos \jjh + 2) - 



2At(l-9)(A<x + Bfi) 



ir- 



' \ — cos \}/h) 



X 4 = 



B-At6(A gl +B g2 ) 



(cosi//A+2) + 



2At6(Aa + B[t) 



1 — cos l///l), 



Al6 = 



2(l+(l-fl)A/g 2 ) 2(\-6)Atp 
3 h 1 



Now on inserting the trial solutions (one Fourier mode out of the 
full solution) at a given point x m O m = A£ J ' exp(imtj/h) and 
jy m = Be; 1 exp(iimj/h) into equation (16) and rearranging the 
equations, where A and B are the harmonics amplitudes, 1// is 
the mode number, h is the element size and i = s/ — 1 we get 



[X 2 + iY\? +1 = [Xi+iY]?, 
[X4 + iY]? +1 = [X 3 + iY]i?, 



(17) 



where 



X V 



A + At(\-e)(Af l+ Bf 2 ) 



(cos \ph + 2) - 



2At(l-9)Aoi 



h 2 



'1 — cosij/h), 



Table 1 1. Maximum error, ratio and order of convergence of 
the proposed scheme for V(x,t). 



f=0.1 






Order of 


f=0.5 




Order of 


n 




Ratio 


Conver. 


U 


Ratio 


Conver. 


16 


1 .26E-02 






4.60E-02 














32 


3.17E-03 


3.9819 


1 .9834 


1.21E-02 


3.7833 


1.9196 


64 


7.95 E-04 


3.9956 


1 .9984 


3.08E-03 


3.9443 


1 .9797 


128 


2.59E-04 


3.9988 


1 .9995 


7.73E-04 


3.9859 


1 .9949 


256 


4.97 E-05 


3.9992 


1 .9997 


1 .93E-04 


3.9961 


1 .9985 


512 


1 .24E-05 


3.9981 


1 .9993 


4.84E-05 


3.9975 


1.9991 


1024 


3.1 1 E-06 


3.9927 


1 .9973 


1.21 E-05 


3.9933 


1 .9975 
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F = 0. 

On direct calculation of equation (17) we obtain 

Xi+iY 
X 2 + iY' 
X 3 + iY 
X 4 + iY' 



(18) 



For stability, the maximum modulus of the eigen-values of the 
matrix has to be less than or equal to one. As 0 = 0.5 is used in the 
proposed scheme, we thus substitute the value of 6 into equation 
(18) and after some algebraic calculation, it can be noticed that 



Icf 



Xf+Y^ 
X 2 +Y 2 ~ ' 
Y 2 



n 



(19) 



X 2 +Y 2 



<1, 



Thus, from (19), the proposed scheme for strongly coupled 
reaction-diffusion equations is unconditionally stable since the 
modulus of the eigen-values must be less than one. This means 
that there are no constraints on grid size h and step size in time 
level At but we should prefer those values ofh and At for which we 
obtain the best accuracy of the scheme. 

Results and Discussion 

To test the accuracy of present method, one example is given in 
this section with L m and relative L 2 error norms are calculated by 



■W 



and in similar way for the numerical solution V(x,t). 

The numerical order of convergence R for both U(x,t) and 
V(x,t) of present method is obtained by using the formula [19]. 
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R = 



Log(L x (ri)/L x (2n)) 
Log(2n/n) 



(20) 



where L^{n) and L aa (2n) are the errors at number of partitions n 
and In respectively. 

We compare the numerical solution obtained by cubic B-spline 
collocation method for strongly coupled reaction-diffusion system 
(1) with known exact solution. 

Example 1 

In this example, we present a strongly coupled reaction- 
diffusion system (1) with numerical solution to show the capability 
and efficiency of cubic B-spline collocation scheme. 

We consider the following problem 



u t = u xx + (2n 2 -\)u-2n 2 v, 0<x< 1, 0< f< T, 



v, = u xx + v xx -v, 
with initial conditions 



0<x<l,0<«<7\ 



(21) 



w(x,0) = sin 2 7rx, 
v(x,0) = cos 2 7ix, 

and boundary conditions as follows: 

f 11,(0,0=0, k»(U)=o, 

\v»(0,f) = 0, v»(l,0 = 0, 



0<x<l, 



0<t<T, 



(22) 



(23) 



It is straightforward to verify that the following are the exact 
solutions [1] 



finite difference scheme in [1]. The numerical errors of the 
proposed method are presented in Table 4 and Table 5 and are 
also depicted graphically in Figure 1 and Figure 2. The absolute 
errors of the numerical solution using finite difference scheme in 
[1] are shown in Table 6 and Table 7 

The solutions are also tabulated in Table 8 and Table 9 with 
different number of partition and different time levels. Results are 
presented graphically for U{x,f) and V(x,i) in Figure 3 and 
Figure 4 for 0 < t < 1 respectively. 

The order of convergence of the present example is calculated 
by the use of the formula given in (20) and which is tabulated in 
Table 10 and Table 11 for U(x,t) and V(x,t) respectively. An 
examination of these tables indicates the method has a nearly 
second order of convergence. 

Concluding Remarks 

In this paper, a numerical method which incorporates a usual 
finite difference scheme with cubic B-spline is presented for solving 
the strongly coupled reaction diffusion system. A finite difference 
approach is used to discretize the time derivatives and cubic B- 
spline is used to interpolate the solutions at each time level. It is 
noted that sometimes the accuracy of solution reduces as time 
increases due to the time truncation errors of time derivative term 
[19]. However the cubic B-spline method used in this work is 
simple and straight forward to apply. The computed results show 
that the cubic B-spline gives reasonable solutions which are 
comparable with finite difference scheme with smaller space steps. 
The obtained solution to the reaction diffusion system for various 
time levels have been compared with the exact solution by finding 
the Lqo and L2 errors. An advantage of using the cubic B-spline 
method outlined in this paper is that it can give accurate solutions 
at any intermediate point in the space direction. 



u(x,t) = e 'sin 2 7ix, 



v(x,t) = e ' cos 2 nx, 



We use the cubic B-spline collocation method (11 ) — ( 1 2) and (14)- 
(15) to compute the numerical solution of (21)— (23). Table 2 and 
Table 3 show the acceptable comparison between numerical and 
exact solutions at some grid points at t = 1 .0 with different number 
of partitions. 

This problem is tested using different values of A and At to show 
the capability of the proposed method for solving the system (21)- 
(23). The final time is taken as T=l. The maximum absolute 
errors of the method at some grid points are comparable with 
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